## libraries
library(DESeq2)
library(limma)
library(edgeR)
library(DT)
library(ComplexHeatmap)
#library(circlize)
library(corrplot)
library(survival)
#library(RColorBrewer)
# gene sets
data(human_gene_signatures)
ind_genes <- human_gene_signatures
# colors
data(color_palettes)
# variables
irf <- "Best Confirmed Overall Response"
ml <- "FMOne mutation burden per MB"
goi <- names(ind_genes)
tablesDir <- system.file("tables",
package="IMvigor210CoreBiologies")
# font sizes
labCex <- 0.9
namesCex <- 0.9
legendCex <- 0.9
titleCex <- 1
axisCex <- 0.9
titleF <- 1
# load
data(cds)
cds2 <- cds
data(fmone)
fmi <- fmone
# normalize
geneNames <- setNames(fData(cds2)$Symbol,
as.character(rownames(fData(cds2))))
voomD <- filterNvoom(counts(cds2),
minSamples=ncol(counts(cds2))/10,
minCpm=0.25)
m <- voomD$E
m <- t(scale( t( m ),
center=TRUE,
scale=TRUE)
)
# add signature scores to pData()
m2 <- m
rownames(m2) <- geneNames[rownames(m2)]
# calculate gene set scores
for (sig in goi) {
pData(cds2)[, sig] <- NA
genes <- ind_genes[[sig]]
genes <- genes[genes %in% rownames(m2)]
tmp <- m2[genes, , drop=FALSE]
pData(cds2)[, sig] <- gsScore(tmp)
}
pData(cds2)$MKI67 <- m2["MKI67", ]
oldMar <- par()$mar
par(mar=c(5.5, 4.1, 2, 0.5))
feat <- "IC Level"
tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]) & pData(cds2)[, irf] != "NE", ]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
ic <- table(tmpDat[, feat], tmpDat[, "binaryResponse"])
pval <- signif(fisher.test(ic)$p.value, 2)
print(paste("Fisher P for IC by binary response:", pval))
## [1] "Fisher P for IC by binary response: 0.0038"
for (rr in c("CR", "PR")) {
tmpDat$group <- ifelse(tmpDat[, irf] == rr, "R", "NR")
tmpDat$group <- factor(tmpDat$group)
ic <- table(tmpDat[, feat], tmpDat$group)
print(rr)
print("all")
print(signif(fisher.test(ic)$p.value, 2))
ic01 <- signif(fisher.test(ic[c("IC0", "IC1"),])$p.value, 2)
ic02 <- signif(fisher.test(ic[c("IC0", "IC2+"),])$p.value, 2)
ic12 <- signif(fisher.test(ic[c("IC1", "IC2+"),])$p.value, 2)
ic2 <- rbind(ic01=colSums(ic[c("IC0", "IC1"),]),
ic2=ic["IC2+",])
ic2 <- signif(fisher.test(ic2)$p.value, 2)
tmp <- c("IC0vsIC1"=ic01, "IC0vs2+"=ic02, "IC1vs2+"=ic12)
tmp <- ifelse(tmp * 3 > 1, 1, tmp * 3)
print(c(tmp, "IC2vsOthers"=ic2))
}
## [1] "CR"
## [1] "all"
## [1] 0.0019
## IC0vsIC1 IC0vs2+ IC1vs2+ IC2vsOthers
## 1.0000 0.0126 0.0171 0.0006
## [1] "PR"
## [1] "all"
## [1] 0.53
## IC0vsIC1 IC0vs2+ IC1vs2+ IC2vsOthers
## 1.00 0.93 1.00 0.30
ic <- table(tmpDat[, feat], tmpDat[, irf])
nSamples <- rowSums(ic)
ic <- prop.table(t(ic),
margin=2)
a <- barplot(ic,
ylab="fraction of patients",
cex.names=namesCex,
cex.axis=axisCex,
cex.lab=labCex,
legend.text=rownames(ic),
col=color_palettes$irf_palette,
width=0.16,
xlim=c(0,1),
args.legend=list(bty="n",
cex=legendCex,
x="topright"),
xaxt="n")
text(x = a,
y = par("usr")[3] - 0.03,
labels = levels(tmpDat[, feat]),
srt = -45,
xpd = TRUE,
adj=0,
cex=namesCex)
mtext(nSamples,
side=3,
at=a,
line=0,
cex=0.7)
mtext(paste("IC", "response", sep=", "),
side=3,
at=a[2],
line=1,
font=titleF,
cex=titleCex)
par(mar=oldMar)
addDots <- FALSE
orgMar <- par()$mar
par(mar=c(5.5, 4.1, 2, 0.5))
tmpDat <- pData(cds2)[pData(cds2)[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
pval <- signif(wilcox.test(tmpDat[, "CD 8 T effector"] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P for Teff by binary response:", pval))
## [1] "Wilcoxon test P for Teff by binary response: 0.0056"
pval <- signif(t.test(tmpDat[, "CD 8 T effector"] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P for Teff by binary response:", pval))
## [1] "T test P for Teff by binary response: 0.0087"
yLab <- "CD8 T-effector score"
a <- boxplot(tmpDat[, "CD 8 T effector"] ~ tmpDat[, irf],
ylab=yLab,
col=color_palettes$irf_palette[levels(tmpDat[, irf])],
cex.main=0.6,
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
ylim=c(-8, 6.5),
whisklty = 1)
if (addDots) {
points(y=tmpDat[, "CD 8 T effector"],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext("CD8 T-effector, response",
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
ind <- tmpDat[, irf] %in% c("CR", "PR")
pvalCRpr <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("CR", "SD")
pvalCRsd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("CR", "PD")
pvalCRpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
pvalsCR <- c(CRpr=pvalCRpr,
CRsd=pvalCRsd,
CRpd=pvalCRpd)
ind <- tmpDat[, irf] %in% c("PR", "SD")
pvalPRsd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("PR", "PD")
pvalPRpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("SD", "PD")
pvalSDpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
pvalsOther <- c(PRsd=pvalPRsd,
PRsd=pvalPRpd,
SDpd=pvalSDpd)
pvalsAll <- c(pvalsCR, pvalsOther)
pvalsAll <- ifelse(pvalsAll * length(pvalsAll) > 1, 1, pvalsAll * length(pvalsAll))
print("t test P for all possible groupings:")
## [1] "t test P for all possible groupings:"
print(pvalsAll)
## CRpr CRsd CRpd PRsd PRsd SDpd
## 0.0386777173 0.0667884675 0.0002873751 1.0000000000 1.0000000000 0.2370005851
pvalsStar <- pLevel(pvalsAll)
segments(x0=1,
x1=4,
y0=6.45,
y1=6.45,
col="black",
xpd=TRUE,
lwd=1)
text(x=2.5,
y=6.75,
labels=pvalsStar["CRpd"],
font=1,
cex=1.2,
xpd=TRUE)
segments(x0=1,
x1=3,
y0=6.06,
y1=6.06,
col="black",
xpd=TRUE,
lwd=1)
text(x=2,
y=6.28,
labels=pvalsStar["CRsd"],
font=1,
cex=1.6,
xpd=TRUE)
segments(x0=1,
x1=2,
y0=5.5,
y1=5.5,
col="black",
xpd=TRUE,
lwd=1)
text(x=1.5,
y=5.75,
labels=pvalsStar["CRpr"],
font=1,
cex=1.2,
xpd=TRUE)
par(mar=orgMar)
crAll <- factor(ifelse(tmpDat[, irf] == "CR", "CR", "nonCR"))
pvalCRall <- t.test(tmpDat[, "CD 8 T effector"] ~ crAll)$p.value
print(paste("t test P for CR vs all:", pvalCRall))
## [1] "t test P for CR vs all: 0.000205031608277682"
ind <- tmpDat[, irf] != "CR"
print("ANOVA P for PR/SD/PD")
## [1] "ANOVA P for PR/SD/PD"
getPfromAnova(tmpDat[ind, "CD 8 T effector"], droplevels(tmpDat[ind, irf]))
## [1] 0.0847
bio <- "CD 8 T effector"
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, bio]),]
tmpDat$group <- cut(tmpDat[, bio],
quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
labels=c("1st Quart",
"2nd Quart",
"3rd Quart",
"4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
plotSurvival2(survFit=fittedSurv,
survDiff=diffSurv,
diff.factor=as.factor(tmpDat$group),
main="CD8 T-effector, survival",
cols=rev(color_palettes[["irf_palette"]][1:4]),
ltypes=rep(1, nlevels(tmpDat$group)),
pval="none",
#coxphData=coxphDat,
plotMedian=TRUE,
plot.nrisk=FALSE,
xLab="OS (months)",
lwd=3,
cexMedian=axisCex - 0.1,
cexLegend=legendCex,
cex.lab=labCex)
print(coxph(Surv(os, censOS) ~
as.integer(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.integer(tmpDat$group) -0.15123 0.85965 0.05803 -2.606 0.00916
##
## Likelihood ratio test=6.81 on 1 df, p=0.009064
## n= 348, number of events= 232
addDots <- FALSE
alternative <- TRUE
## response by Mutation burden
oldMar <- par()$mar
par(mar=c(5.5, 4.1, 2, 0.5))
bio <- ml
feat <- irf
clinDat2 <- pData(cds2)[!is.na(pData(cds2)[, bio]) & pData(cds2)[, bio] > 0,]
tmpDat <- clinDat2[!is.na(clinDat2[, feat]),]
tmpDat <- tmpDat[tmpDat[, feat] != "NE",]
tmpDat[, feat] <- droplevels(tmpDat[, feat])
ind <- tmpDat[, bio] > 0
tmpDat[, bio] <- log2(tmpDat[, bio])
yLab <- "TMB"
a <- boxplot(tmpDat[ind, bio] ~ tmpDat[ind, feat],
ylab=yLab,
col=color_palettes$irf_palette,
cex.main=0.6,
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
whisklty = 1,
ylim=c(1,log2(80)),
yaxt="n")
axis(2,
at=log2(c(1,2,5,10,20,50)),
labels=c(1,2,5,10,20,50),
cex=axisCex)
if (addDots) {
points(y=tmpDat[ind, bio],
x=jitter(as.numeric(tmpDat[ind, feat]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, feat]),
line=0,
cex=0.8)
pval <- signif(wilcox.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon Test p TMB by response:", pval))
## [1] "Wilcoxon Test p TMB by response: 1.4e-07"
pval <- signif(t.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value,
2)
print(paste("T Test p TMB by response:", pval))
## [1] "T Test p TMB by response: 6.9e-07"
mtext("TMB, response",
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
if (alternative) {
yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
x1=2,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3,
x1=4,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=1.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3.5,
x1=3.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=3.5,
y0=par("usr")[4] - yunit * 2,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
text(y=par("usr")[4] - yunit ,
x=2.5,
labels=pLevel(pval),
cex=axisCex)
}
par(mar=oldMar)
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, ml]),]
tmpDat$group <- cut(tmpDat[, ml],
quantile(tmpDat[, ml], c(0, 0.25, 0.5, 0.75, 1)),
include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
labels=c("1st Quart",
"2nd Quart",
"3rd Quart",
"4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
plotSurvival2(survFit=fittedSurv,
survDiff=diffSurv,
diff.factor=as.factor(tmpDat$group),
main="TMB, survival",
cols=rev(color_palettes[["irf_palette"]][1:4]),
ltypes=rep(1, nlevels(tmpDat$group)),
pval="none",
#coxphData=coxphDat,
plotMedian=TRUE,
plot.nrisk=FALSE,
xLab="OS (months)",
lwd=3,
cexMedian=axisCex,
cexLegend=legendCex,
cex.lab=labCex)
print(coxph(Surv(os, censOS) ~
as.integer(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.integer(tmpDat$group) -0.29531 0.74430 0.06924 -4.265 2e-05
##
## Likelihood ratio test=19 on 1 df, p=1.309e-05
## n= 272, number of events= 173
ind <- !is.na(pData(fmi)$binaryResponse)
pData(fmi)$binaryResponse <- factor(pData(fmi)$binaryResponse,
levels=c("CR/PR", "SD/PD"))
par(mar=c(5.5, 4.1, 2, 0))
path <- "DDR"
ids <- ind_genes[[path]]
ids <- ids[!ids %in% "TP53"]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)
pval <- signif(fisher.test(table(pData(fmi)[ind, "binaryResponse"], mutStatus[ind] > 0))$p.value)
print(paste("Fisher P binary response by DDR mutation status:", pval))
## [1] "Fisher P binary response by DDR mutation status: 0.0116883"
b <- table(responder = fmi$binaryResponse[ind],
mutant = factor(mutStatus[ind] > 0))
nSamples <- colSums(b)
d <- prop.table(b,
margin=2)
a <- barplot(d,
legend.text=colnames(d),
col=color_palettes$response_palette[rownames(d)],
cex.names=namesCex,
cex.axis=axisCex,
cex.lab=labCex,
width=0.2,
xlim=c(0,1),
ylab="fraction of patients",
xaxt="n",
args.legend=list(x="topright",
bty="n",
fill=color_palettes$response_palette,
legend=c("CR/PR", "SD/PD"),
cex=legendCex))
text(x = a,
y = par("usr")[3] - 0.06,
labels = ifelse(colnames(d) == FALSE, "non-mutant", "mutant"),
srt = -45,
xpd = TRUE,
adj=0,
cex=namesCex,
xpd=TRUE)
mtext(nSamples,
side=3,
at=a,
line=0,
cex=0.7)
mtext(paste(path, "w/o TP53, response"),
side=3,
at=mean(a),
line=1,
font=titleF,
cex=titleCex)
mtext(pLevel(pval),
at=mean(a),
side=1,
line=0.05,
cex=namesCex,
xpd=TRUE)
segments(x0=a[1] - 0.05,
x1=a[2] + 0.05,
y0=-0.02,
y1=-0.02,
col="black",
xpd=TRUE,
lwd=1)
par(mar=oldMar)
ind <- !is.na(pData(fmi)$binaryResponse)
pathFMIout <- lapply(c("DDR", "CC Reg"), function(path) {
ids <- ind_genes[[path]]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)
t <- table(responder = fmi$binaryResponse[ind],
mutant = factor(mutStatus[ind] > 0))
f <- fisher.test( t )
setNames(c( path, sum(mutStatus[ind] > 0), "response", f$estimate, f$p.value ),
c("Pathway", "n mutant", "category", "estimate", "PVal"))
})
pathFMIout <- do.call(rbind, pathFMIout)
pathFMIout <- as.data.frame(pathFMIout,
stringsAsFactors=FALSE)
pathFMIout$estimate <- signif(as.numeric(pathFMIout$estimate))
pathFMIout$PVal <- signif(as.numeric(pathFMIout$PVal))
#pathFMIout$AdjP <- p.adjust(pathFMIout$PVal, method="BH")
ind <- !is.na(pData(fmi)[,ml])
pathFMIout2 <- lapply(c("DDR", "CC Reg"), function(path) {
ids <- ind_genes[[path]]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)
f <- t.test( as.numeric(pData(fmi)[ind,ml]) ~ mutStatus[ind] > 0,
conf.int = TRUE )
setNames(c( path, sum(mutStatus[ind] > 0), "TMB", f$estimate[1] - f$estimate[2], f$p.value ),
c("Pathway", "n mutant", "category", "estimate", "PVal"))
})
pathFMIout2 <- do.call(rbind, pathFMIout2)
pathFMIout2 <- as.data.frame(pathFMIout2,
stringsAsFactors=FALSE)
pathFMIout2$estimate <- signif(as.numeric(pathFMIout2$estimate))
pathFMIout2$PVal <- signif(as.numeric(pathFMIout2$PVal))
#pathFMIout2$AdjP <- p.adjust(pathFMIout2$PVal, method="BH")
pathFMIout <- rbind(pathFMIout, pathFMIout2)
datatable(pathFMIout,
rownames=FALSE)
tmp <- any_mutation(fmi)
tmp <- tmp[rowSums(tmp) > 0, ]
pathFMIout3 <- lapply(rownames(tmp), function(gene) {
mutStatus <- tmp[gene,]
ind <- !is.na(pData(fmi)[,ml])
f <- try(t.test( as.numeric(pData(fmi)[ind,ml]) ~ mutStatus[ind] > 0,
conf.int = TRUE ))
if (!is(f, "try-error")) {
a <- setNames(c( gene,
sum(mutStatus[ind] > 0),
"TMB",
f$estimate[1]-f$estimate[2],
f$p.value ),
c("Gene", "n mutant", "category", "estimate", "PVal"))
} else {
a <- setNames(c( gene,
sum(mutStatus[ind] > 0),
"TMB",
NA,
NA ),
c("Gene", "n mutant", "category", "estimate", "PVal"))
}
ind <- !is.na(pData(fmi)$binaryResponse)
t <- table(responder = fmi$binaryResponse[ind],
mutant = factor(mutStatus[ind] > 0, levels=c("TRUE", "FALSE")))
f <- fisher.test( t )
b <- setNames(c( gene, sum(mutStatus[ind] > 0), "response", f$estimate, f$p.value ),
c("Gene", "n mutant", "category", "estimate", "PVal"))
rbind(a, b)
})
Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations
pathFMIout3 <- do.call(rbind, pathFMIout3)
pathFMIout3 <- as.data.frame(pathFMIout3,
stringsAsFactors=FALSE)
pathFMIout3$estimate <- signif(as.numeric(pathFMIout3$estimate))
pathFMIout3$PVal <- signif(as.numeric(pathFMIout3$PVal))
## adjust p-values
pathFMIout3$AdjP <- NA
pathFMIout3$AdjP[pathFMIout3$category == "response"] <- p.adjust(
pathFMIout3$PVal[pathFMIout3$category == "response"] , method="BH")
pathFMIout3$AdjP[pathFMIout3$category == "TMB"] <- p.adjust(
pathFMIout3$PVal[pathFMIout3$category == "TMB"] , method="BH")
pathFMIout3$inDDR <- pathFMIout3$Gene %in% ind_genes$DDR
pathFMIout3$inCCReg <- pathFMIout3$Gene %in% ind_genes$"CC Reg"
datatable(pathFMIout3,
rownames=FALSE)
path <- "DDR"
fmiPathP <- pathFMIout[pathFMIout$category == "TMB", ]
fmiGeneP <- pathFMIout3[pathFMIout3$category == "TMB", ]
ids <- ind_genes[[path]]
ids <- ids[ids %in% featureNames(fmi)]
fmiGoi <- fmi[featureNames(fmi) %in% ids, ]
sOrder <- order(as.numeric(pData(fmiGoi)[, ml]), decreasing=TRUE)
mutGoi <- any_mutation(fmiGoi)
tmp1 <- mutGoi[rowSums(mutGoi) > 0,]
tmp2 <- colSums(tmp1)
tmp2 <- tmp2 > 0
if ("TP53" %in% ids) {
tmp3 <- mutGoi[ids[! ids %in% "TP53"], ]
tmp3 <- colSums(tmp3)
tmp3 <- tmp3 > 0
tmp2 <- rbind(tmp2, tmp3)
}
toSort <- rowSums(tmp1)
toSort <- sort(toSort, decreasing=TRUE)
mutGoi <- rbind(tmp1[names(toSort),], tmp2)
rownames(mutGoi)[rownames(mutGoi) == "tmp2"] <- "pathway"
if ("TP53" %in% ids) {
rownames(mutGoi)[rownames(mutGoi) == "tmp3"] <- "pathway w/o TP53"
}
mutGoi <- apply(mutGoi, 2, function(x) {
ifelse(x, "mutant", "")
})
fmiCols <- c(mutant="black")
alter_fun = function(x, y, w, h, v) {
if(v["mutant"]) grid.rect(x, y, w*0.9, h*0.9, gp = gpar(fill = fmiCols["mutant"],
col = NA))
}
ha = HeatmapAnnotation(TMB = anno_barplot(as.numeric(pData(fmiGoi)[sOrder, ml]),
border=FALSE,
gp = gpar(fill="black")),
annotation_legend_param=list(gp = gpar(fontsize = 1)),
show_annotation_name = TRUE,
annotation_name_side="left",
annotation_name_gp = gpar(fontsize = 10))
rownames(mutGoi) <- paste0(rownames(mutGoi),
ifelse(rownames(mutGoi) %in% fmiGeneP$Gene[fmiGeneP$AdjP < 0.05], "*", ""))
op1 <- oncoPrint(mutGoi[, sOrder],
col=fmiCols,
alter_fun=alter_fun,
row_order=NULL,
column_order=NULL,
column_title = "Mutations in cell cycle regulation genes",
top_annotation=ha,
show_pct=TRUE,
pct_gp = gpar(fontsize = 9),
row_names_gp = gpar(fontsize = 9),
column_title_gp = gpar(fontsize = 10),
heatmap_legend_param=list(title ="Mutation Status"),
show_row_barplot = FALSE,
show_heatmap_legend = FALSE)
op2 <- oncoPrint(mutGoi[, sOrder],
col=fmiCols,
alter_fun=alter_fun,
row_order=NULL,
column_order=NULL,
column_title = "Mutations in DDR genes",
top_annotation=ha,
show_pct=TRUE,
pct_gp = gpar(fontsize = 9),
row_names_gp = gpar(fontsize = 9),
column_title_gp = gpar(fontsize = 10),
heatmap_legend_param=list(title ="Mutation Status"),
show_row_barplot = FALSE,
show_heatmap_legend = FALSE)
draw(op2,
padding = unit(c(1, 5, 0, 1), "mm"),
row_title="mutation prevalence",
row_title_gp=gpar(fontsize = 10))
plotDots <- FALSE
binary <- FALSE
orgMar <- par()$mar
tmpDat <- pData(cds2)
tmpDat$TGFB1 <- scale(m2["TGFB1", ], center=TRUE, scale=TRUE)
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
sig <- "TGFB1"
par(mar=c(5.5, 4.1, 2, 0.5))
pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))
## [1] "Wilcoxon test P TGFB1 score by binary Response: 5.6e-05"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P", sig, "score by binary Response:", pval))
## [1] "T test P TGFB1 score by binary Response: 0.00011"
yLab <- paste(sig, "expression")
a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
ylab=yLab,
col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
ylim=c(-5, 5),
whisklty = 1)
if (plotDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext(paste0(sig, ", response"),
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
x1=2,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3,
x1=4,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=1.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3.5,
x1=3.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=3.5,
y0=par("usr")[4] - yunit * 2,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
text(y=par("usr")[4] - yunit ,
x=2.5,
labels=pLevel(pval),
cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex))
par(mar=orgMar)
bio <- "TGFB1"
tmpDat <- pData(cds2)
tmpDat$TGFB1 <- m2["TGFB1", ]
tmpDat <- tmpDat[!is.na(tmpDat[, bio]),]
tmpDat$group <- cut(tmpDat[, bio],
quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
labels=c("1st Quart",
"2nd Quart",
"3rd Quart",
"4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
plotSurvival2(survFit=fittedSurv,
survDiff=diffSurv,
diff.factor=as.factor(tmpDat$group),
main="TGFB1, survival",
cols=rev(color_palettes[["irf_palette"]][1:4]),
ltypes=rep(1, nlevels(tmpDat$group)),
pval="none",
#coxphData=coxphDat,
plotMedian=TRUE,
plot.nrisk=FALSE,
xLab="OS (months)",
lwd=3,
cexMedian=axisCex - 0.1,
cexLegend=legendCex,
cex.lab=labCex)
print(coxph(Surv(os, censOS) ~
as.integer(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.integer(tmpDat$group) 0.15298 1.16530 0.05905 2.591 0.00958
##
## Likelihood ratio test=6.74 on 1 df, p=0.009433
## n= 348, number of events= 232
table(IC=pData(cds2)[, "IC Level"], TC=pData(cds2)[, "TC Level"])
## TC
## IC TC0 TC1 TC2+
## IC0 93 4 0
## IC1 111 10 11
## IC2+ 71 8 39
tmp <- pData(cds2)[!is.na(pData(cds2)[, "TC Level"]), ]
a <- ifelse(tmp[, "TC Level"] %in% "TC0", "TC-", "TC+")
b <- ifelse(tmp[, "IC Level"] %in% "IC0", "IC-", "IC+")
table(IC=b, TC=a)
## TC
## IC TC- TC+
## IC- 93 4
## IC+ 182 68
prop.table(table(IC=b, TC=a), 2)
## TC
## IC TC- TC+
## IC- 0.33818182 0.05555556
## IC+ 0.66181818 0.94444444
par(mar=c(5.5, 4.1, 2, 0.5))
for (feat in "TC Level") {
tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]) & pData(cds2)[, irf] != "NE", ]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
ic <- table(tmpDat[, feat], tmpDat[, "binaryResponse"])
pval <- signif(fisher.test(ic)$p.value, 2)
print(paste("Fisher P for TC by binary response:", pval))
ic <- table(tmpDat[, feat], tmpDat[, irf])
nSamples <- rowSums(ic)
ic <- prop.table(t(ic),
margin=2)
a <- barplot(ic,
#main=paste("P:", pval),
ylab="fraction of patients",
#xlab=feat,
cex.names=namesCex,
cex.axis=axisCex,
cex.lab=labCex,
legend.text=rownames(ic),
col=color_palettes[["irf_palette"]],
width=0.16,
xlim=c(0,1),
#las=2,
args.legend=list(bty="n",
#bg=alpha("white", 0.7),
#box.col=alpha("white", 0.7),
cex=1.1,
x="topright"),
xaxt="n")
text(x = a,
y = par("usr")[3] - 0.03,
labels = levels(tmpDat[, feat]),
srt = -45,
#pos = 1,
xpd = TRUE,
adj=0,
cex=namesCex)
mtext(nSamples,
side=3,
at=a,
line=0,
cex=0.7)
#mtext(paste("P:", pval),
# side=1,
# at=a[2],
# line=3.5,
# cex=0.8)
mtext(paste("TC", "response", sep=", "),
side=3,
at=a[2],
line=1,
font=titleF,
cex=titleCex)
}
## [1] "Fisher P for TC by binary response: 0.72"
oldMar <- par()$mar
par(mar=c(5.5, 3.5, 2, 0.5))
icAss <- read.csv(file.path(tablesDir,
"SupplementaryTableS10.csv"),
stringsAsFactor=FALSE,
check.names=FALSE)
icAss <- icAss[, 1:8]
icAss <- icAss[!is.na(icAss$Symbol), ]
interferome <- read.table(file.path(tablesDir,
"20170503081937_GeneSearchResults_InterferomeTop300IC.txt"),
stringsAsFactor=FALSE,
sep="\t",
check.names=FALSE,
skip=19,
header=TRUE)
cols <- "darkgrey"
plot(x=icAss[,"logFC"],
y=-log10(icAss[,"P (adj.)"]),
pch=16,
cex=0.8,
col=cols,
ylab="",
xlab="",
lwd=2,
cex.axis=axisCex,
cex.lab=labCex)
mtext(expression(paste("-log"["10"], " adj. p")),
side=2,
line=2,
cex=axisCex)
abline(h=-log10(0.05), lty=2, col="black")
gs <- list("Interferome"=interferome$"Gene Name"[interferome$"Gene Name" %in% icAss$Symbol[1:307]],
"CD 8 T effector"=ind_genes[["CD 8 T effector"]],
"Immune Checkpoint"=ind_genes[["Immune Checkpoint"]])
gsC <- list("Interferome"="orange1",
"CD 8 T effector"="magenta2",
"Immune Checkpoint"="olivedrab3")
for (g in names(gs)) {
ind <- which(icAss$Symbol %in% gs[[g]])
points(x=icAss[ind,"logFC"],
y=-log10(icAss[ind,"P (adj.)"]),
pch=16,
cex=0.8,
col=gsC[[g]])
}
legend("topleft",
pch=16,
col=c("orange1", "magenta2", "olivedrab3"),
legend=c("IFNg inducible", "CD8 T-effector", "Immune checkpoint"),
bty = "n",
cex=0.6)
labelSub <- icAss[icAss$Symbol %in%
c(ind_genes[["CD 8 T effector"]], ind_genes[["Immune Checkpoint"]]), ]
labelPos <- c(2,
2,
2,
4,
rep(3, 2),
4,
3,
4,
2,
1,
2,
1,
2,
1)
text(x=labelSub[,"logFC"],
y=-log10(labelSub[,"P (adj.)"]),
labels=labelSub$Symbol,
cex= 0.4,
pos=labelPos)
mtext("Expression, IC",
side=3,
at=0.5,
line=1,
font=titleF,
cex=titleCex)
mtext("log fold change",
side=1,
at=0.5,
line=2,
cex=axisCex)
par(mar=oldMar)
addDots <- FALSE
sig <- "CD 8 T effector"
feat <- "IC Level"
oldMar <- par()$mar
par(mar=c(5.5, 4.1, 2, 2))
tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]), ]
a <- boxplot(tmpDat[, sig] ~ tmpDat[, feat],
ylab="CD8 T-effector score",
col=color_palettes$ic_palette[levels(tmpDat[, feat])],
cex.axis=axisCex,
cex.lab=labCex,
cex.names=namesCex,
whisklty = 1,
xaxt="n")
axis(1,
at=1:nlevels(tmpDat[, feat]),
labels=rep("", nlevels(tmpDat[, feat])))
text(x = 1:nlevels(tmpDat[, feat]),
y = par("usr")[3] - 1.5,
labels = levels(tmpDat[, feat]),
xpd = TRUE,
cex=namesCex)
if (addDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, feat]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, feat]),
line=0,
cex=0.8)
mtext(paste("CD8 T-effector", "IC", sep=", "),
side=3,
at=2,
line=1,
font=titleF,
cex=titleCex)
pval <- signif(getPfromAnova(tmpDat[, sig], tmpDat[, feat]), 2)
print(paste("ANOVA P for Teff signature by IC level:", pval))
## [1] "ANOVA P for Teff signature by IC level: 4.2e-35"
par(mar=oldMar)
addDots <- FALSE
alternative <- TRUE
## response by Mutation burden
oldMar <- par()$mar
par(mar=c(5.5, 4.1, 2, 0.5))
bio <- "Neoantigen burden per MB"
feat <- irf
clinDat2 <- pData(cds2)[!is.na(pData(cds2)[, bio]) & pData(cds2)[, bio] > 0,]
tmpDat <- clinDat2[!is.na(clinDat2[, feat]),]
tmpDat <- tmpDat[tmpDat[, feat] != "NE",]
tmpDat[, feat] <- droplevels(tmpDat[, feat])
ind <- tmpDat[, bio] > 0
tmpDat[, bio] <- log2(tmpDat[, bio])
yLab <- "TNB"
a <- boxplot(tmpDat[ind, bio] ~ tmpDat[ind, feat],
ylab=yLab,
col=color_palettes$irf_palette,
cex.main=0.6,
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
whisklty = 1,
ylim=c(-5, log2(30)),
yaxt="n")
axis(2,
at=log2(c(0.05,0.20,1.00,5.00,20)),
labels=round(c(0.05,0.20,1.00,5.00,20.00), 2),
cex=axisCex)
if (addDots) {
points(y=tmpDat[ind, bio],
x=jitter(as.numeric(tmpDat[ind, feat]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, feat]),
line=0,
cex=0.8)
pval <- signif(wilcox.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon p-value, neoantigens by response:", pval))
## [1] "Wilcoxon p-value, neoantigens by response: 5.3e-09"
pval <- signif(t.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value,
2)
print(paste("T test p-value, neoantigens by response:", pval))
## [1] "T test p-value, neoantigens by response: 2.7e-09"
mtext("Neoantigens, response",
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
if (alternative) {
yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
x1=2,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3,
x1=4,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=1.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3.5,
x1=3.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=3.5,
y0=par("usr")[4] - yunit * 2,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
text(y=par("usr")[4] - yunit ,
x=2.5,
labels=pLevel(pval),
cex=axisCex)
}
par(mar=oldMar)
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, "Neoantigen burden per MB"]),]
tmpDat$group <- cut(tmpDat[, "Neoantigen burden per MB"],
quantile(tmpDat[, "Neoantigen burden per MB"], c(0, 0.25, 0.5, 0.75, 1)),
include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
labels=c("1st Quart",
"2nd Quart",
"3rd Quart",
"4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
plotSurvival2(survFit=fittedSurv,
survDiff=diffSurv,
diff.factor=as.factor(tmpDat$group),
main="TNB, survival",
cols=rev(color_palettes[["irf_palette"]][1:4]),
ltypes=rep(1, nlevels(tmpDat$group)),
pval="none",
#coxphData=coxphDat,
plotMedian=TRUE,
plot.nrisk=FALSE,
xLab="OS (months)",
lwd=3,
cexMedian=axisCex - 0.1,
cexLegend=legendCex,
cex.lab=labCex)
print(coxph(Surv(os, censOS) ~
as.integer(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.integer(tmpDat$group) -0.26824 0.76472 0.07072 -3.793 0.000149
##
## Likelihood ratio test=14.75 on 1 df, p=0.0001229
## n= 245, number of events= 160
print(coxph(Surv(os, censOS) ~
as.ordered(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.ordered(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.ordered(tmpDat$group).L -0.6645 0.5145 0.1735 -3.829 0.000129
## as.ordered(tmpDat$group).Q -0.3235 0.7236 0.1649 -1.962 0.049800
## as.ordered(tmpDat$group).C -0.2942 0.7451 0.1566 -1.878 0.060337
##
## Likelihood ratio test=21.28 on 3 df, p=9.203e-05
## n= 245, number of events= 160
keggRes <- read.csv(file.path(tablesDir,
"SupplementaryTableS1.csv"),
stringsAsFactor=FALSE,
check.names=FALSE)
keggRes <- keggRes[, 1:6]
keggRes <- keggRes[!is.na(keggRes$S), ]
ind <- rev(1:nrow(keggRes))
keggRes <- keggRes[ind,]
up <- keggRes$Direction == "Up"
keggRes$"P (adj.)"[up] <- log10(keggRes$"P (adj.)"[up]) * -1
down <- keggRes$Direction == "Down"
keggRes$"P (adj.)"[down] <- log10(keggRes$"P (adj.)"[down])
resC1 <- read.csv(file.path(tablesDir,
"SupplementaryTableS2.csv"),
stringsAsFactor=FALSE,
check.names=FALSE)
keggRes$RepGenes <- NA_character_
for (i in 1:nrow(keggRes)) {
intGenes <- unlist(strsplit(keggRes$"Identified Genes"[i], ", "))
if (length(intGenes) > 7) {
intGenes <- resC1$Symbol[resC1$Symbol %in% intGenes][1:7]
intGenes <- c(intGenes, "...")
}
intGenes <- paste(intGenes, collapse=", ")
keggRes$RepGenes[i] <- intGenes
}
layout(matrix(c(1,1,0,0,0,0,1,1,0,0,0,0), 2, 6, byrow=TRUE))
proliferationPathways <- c("DNA replication",
"Cell cycle",
"Systemic lupus erythematosus",
"Alcoholism",
"Oocyte meiosis",
"Viral carcinogenesis",
"Progesterone-mediated oocyte maturation")
ddrPathways <- c("Fanconi anemia pathway", "Homologous recombination", "Mismatch repair", "Nucleotide excision repair", "Base excision repair", "p53 signaling pathway")
setCols <- ifelse(keggRes$Name %in% proliferationPathways,
"turquoise4",
ifelse(keggRes$Name %in% ddrPathways,
"magenta4",
ifelse(keggRes$Name %in% "Pathways in cancer", "orange2", "darkgrey")))
a <- barplot(keggRes$"P (adj.)",
horiz=TRUE,
xlab="",
xlim=c(-2, 6),
col=setCols,
axes=FALSE,
cex.main=0.9)
mtext("Kegg pathway",
side=4,
at=max(a[,1]) + 1.5,
las=2,
line=3,
cex=titleCex)
mtext("Top-scoring genes",
side=4,
at=max(a[,1]) + 1.5,
las=2,
line=24,
cex=titleCex)
mtext(side=4,
at=a[,1],
text=paste0(keggRes$Name,
ifelse(keggRes$rep, "*", "")),
las=2,
line=1.5,
cex=(namesCex - 0.1))
mtext(keggRes$RepGenes,
side=4,
at=a[,1],
line=22,
cex=(namesCex - 0.1),
las=2)
axis(side=1,
at=c(-2, seq(from=2, to=6, by=2)),
labels=c(2, seq(from=2, to=6, by=2)),
tick=TRUE,
xpd=TRUE,
cex.axis=0.9)
segments(x0=0,
x1=0,
y0=-0.6,
y1=max(a[,1])+0.7,
col="black",
xpd=TRUE,
lwd=1)
mtext(expression(paste("-log"["10"], " adj. p")),
at=2,
side=1,
line=2.4,
cex=0.9)
mtext(c("Positive association with TMB", "Negative"),
side=2,
line=1.2,
at=c(a[,1][9], mean(c(a[,1][1], a[,1][2])) + 0.2),
cex=axisCex)
segments(x0=-2.1,
x1=-2.1,
y0=a[,1][4] - 0.3,
y1=a[,1][length(a[,1])] + 0.3,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=-2.1,
x1=-2.1,
y0=a[,1][1] - 0.5,
y1=a[,1][3] + 0.5,
col="black",
xpd=TRUE,
lwd=1)
legend("bottomright",
fill=c("turquoise4", "magenta4", "orange2"),
legend=c("Proliferation", "DNA damage response", expression(paste("TGF-", beta))),
bty="n",
cex=legendCex + 0.3)
goi <- c("CD 8 T effector",
"gene19",
"Fanconi",
"MKI67",
"Cell cycle",
"DNA replication",
"Histones",
"Homologous recombination",
"Mismatch Repair",
"Nucleotide excision repair",
"CC Reg")
pLabels <- sub("CC Reg", "Cell cycle regulation", goi)
pLabels <- sub("CD 8 T e", "CD8 T-e", pLabels)
tmp <- pData(cds2)[, goi]
colnames(tmp) <- sub("CC Reg", "Cell cycle regulation", colnames(tmp))
colnames(tmp) <- sub("CD 8 T e", "CD8 T-e", colnames(tmp))
colnames(tmp) <- sub("CC Reg", "Cell cycle regulation", colnames(tmp))
colnames(tmp) <- sub("Fanconi", "Fanconi anemia pathway", colnames(tmp))
colnames(tmp) <- sub("p53", "p53-like", colnames(tmp))
colnames(tmp) <- sub("gene19", "Pan-F-TBRS", colnames(tmp))
M <- cor(tmp)
corrplot.mixed(M,
lower="circle",
upper="number",
number.cex = 0.4,
order="AOE",
tl.cex = 0.9,
tl.col = "black",
xpd=TRUE,
tl.pos="lt",
diag="n")
addDots <- FALSE
sig <- "APOBEC3A"
oldMar <- par()$mar
## plot APOBEC expression by response
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P APOBEC3A expression by binary response:", pval))
## [1] "Wilcoxon test P APOBEC3A expression by binary response: 0.0071"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P APOBEC3A expression by binary response:", pval))
## [1] "T test P APOBEC3A expression by binary response: 0.015"
yLab <- expression(paste("APOBEC3A expression (log"["2"], ")"))
par(mar=c(5.5, 4, 2, 0.5))
a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
ylab=yLab,
col=color_palettes$irf_palette[levels(tmpDat[, irf])],
cex.main=0.6,
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
whisklty = 1)
if (addDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext(paste(sig, "response", sep=", "),
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
## plot APOBEC expression by TMB
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
ind <- tmpDat[, ml] > 0
tmpDat[, ml] <- log2(tmpDat[, ml])
plot(tmpDat[ind, sig] ~ tmpDat[ind, ml],
xlab="",
ylab=expression(paste("APOBEC3A expression (log"["2"], ")")),
pch=16,
col="darkgrey",
cex.axis=0.8,
cex.lab=labCex,
xaxt="n")
axis(1,
at=log2(c(1,2,5,10,20,50)),
labels=c(1,2,5,10,20,50),
cex=axisCex)
fit <- lm(tmpDat[ind, sig] ~ tmpDat[ind, ml])
abline(fit, lty=2)
tmp <- cor.test(tmpDat[ind, sig], tmpDat[ind, ml], method="pearson")
text(paste("cor:", round(tmp$estimate, 2),
"\nP:", signif(tmp$p.value, 2)),
font=2,
cex=0.7,
x=log2(2),
y=3.2)
mtext("TMB",
at=3,
side=1,
line=2.5,
cex=axisCex)
mtext(paste("TMB", sig, sep=", "),
at=3,
side=3,
line=1,
cex=titleCex,
font=titleF)
par(mar=oldMar)
addDots <- FALSE
sig <- "APOBEC3B"
oldMar <- par()$mar
## plot APOBEC expression by response
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P APOBEC3B expression by binary response:", pval))
## [1] "Wilcoxon test P APOBEC3B expression by binary response: 0.0066"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P APOBEC3B expression by binary response:", pval))
## [1] "T test P APOBEC3B expression by binary response: 0.0025"
yLab <- expression(paste("APOBEC3B expression (log"["2"], ")"))
par(mar=c(5.5, 4, 2, 0.5))
a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
ylab=yLab,
col=color_palettes$irf_palette[levels(tmpDat[, irf])],
cex.main=0.6,
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
whisklty = 1)
if (addDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext(paste(sig, "response", sep=", "),
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
## plot APOBEC expression by TMB
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
ind <- tmpDat[, ml] > 0
tmpDat[, ml] <- log2(tmpDat[, ml])
plot(tmpDat[ind, sig] ~ tmpDat[ind, ml],
xlab="",
ylab=expression(paste("APOBEC3B expression (log"["2"], ")")),
pch=16,
col="darkgrey",
cex.axis=0.8,
cex.lab=labCex,
xaxt="n")
axis(1,
at=log2(c(1,2,5,10,20,50)),
labels=c(1,2,5,10,20,50),
cex=axisCex)
fit <- lm(tmpDat[ind, sig] ~ tmpDat[ind, ml])
abline(fit, lty=2)
tmp <- cor.test(tmpDat[ind, sig], tmpDat[ind, ml], method="pearson")
text(paste("cor:", round(tmp$estimate, 2),
"\nP:", signif(tmp$p.value, 2)),
font=2,
cex=0.7,
x=log2(2),
y=2)
mtext("TMB",
at=3,
side=1,
line=2.5,
cex=axisCex)
mtext(paste("TMB", sig, sep=", "),
at=3,
side=3,
line=1,
cex=titleCex,
font=titleF)
par(mar=oldMar)
path <- "CC Reg"
fmiPathP <- pathFMIout[pathFMIout$category == "TMB", ]
fmiGeneP <- pathFMIout3[pathFMIout3$category == "TMB", ]
ids <- ind_genes[[path]]
ids <- ids[ids %in% featureNames(fmi)]
fmiGoi <- fmi[featureNames(fmi) %in% ids, ]
sOrder <- order(as.numeric(pData(fmiGoi)[, ml]), decreasing=TRUE)
mutGoi <- any_mutation(fmiGoi)
tmp1 <- mutGoi[rowSums(mutGoi) > 0,]
tmp2 <- colSums(tmp1)
tmp2 <- tmp2 > 0
if ("TP53" %in% ids) {
tmp3 <- mutGoi[ids[! ids %in% "TP53"], ]
tmp3 <- colSums(tmp3)
tmp3 <- tmp3 > 0
tmp2 <- rbind(tmp2, tmp3)
}
toSort <- rowSums(tmp1)
toSort <- sort(toSort, decreasing=TRUE)
mutGoi <- rbind(tmp1[names(toSort),], tmp2)
rownames(mutGoi)[rownames(mutGoi) == "tmp2"] <- "pathway"
if ("TP53" %in% ids) {
rownames(mutGoi)[rownames(mutGoi) == "tmp3"] <- "pathway w/o TP53"
}
mutGoi <- apply(mutGoi, 2, function(x) {
ifelse(x, "mutant", "")
})
fmiCols <- c(mutant="black")
alter_fun = function(x, y, w, h, v) {
if(v["mutant"]) grid.rect(x, y, w*0.9, h*0.9, gp = gpar(fill = fmiCols["mutant"],
col = NA))
}
ha = HeatmapAnnotation(TMB = anno_barplot(as.numeric(pData(fmiGoi)[sOrder, ml]),
border=FALSE,
gp = gpar(fill="black")),
annotation_legend_param=list(gp = gpar(fontsize = 1)),
show_annotation_name = TRUE,
annotation_name_side="left",
annotation_name_gp = gpar(fontsize = 10))
rownames(mutGoi) <- paste0(rownames(mutGoi),
ifelse(rownames(mutGoi) %in% fmiGeneP$Gene[fmiGeneP$AdjP < 0.05], "*", ""))
op1 <- oncoPrint(mutGoi[, sOrder],
col=fmiCols,
alter_fun=alter_fun,
row_order=NULL,
column_order=NULL,
column_title = "Mutations in cell cycle regulation genes",
top_annotation=ha,
show_pct=TRUE,
pct_gp = gpar(fontsize = 9),
row_names_gp = gpar(fontsize = 9),
column_title_gp = gpar(fontsize = 10),
heatmap_legend_param=list(title ="Mutation Status"),
show_row_barplot = FALSE,
show_heatmap_legend = FALSE)
draw(op1,
padding = unit(c(1, 5, 0, 1), "mm"),
row_title="mutation prevalence",
row_title_gp=gpar(fontsize = 10))
ind <- !is.na(pData(fmi)$binaryResponse)
pData(fmi)$binaryResponse <- factor(pData(fmi)$binaryResponse,
levels=c("CR/PR", "SD/PD"))
par(mar=c(5.5, 4.1, 2, 0))
path <- "CC Reg"
ids <- ind_genes[[path]]
ids <- ids[!ids %in% "TP53"]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)
pval <- signif(fisher.test(table(pData(fmi)[ind, "binaryResponse"], mutStatus[ind] > 0))$p.value)
print(paste("Fisher P binary response by CCR mutation status:", pval))
## [1] "Fisher P binary response by CCR mutation status: 0.31104"
b <- table(responder = fmi$binaryResponse[ind],
mutant = factor(mutStatus[ind] > 0))
nSamples <- colSums(b)
d <- prop.table(b,
margin=2)
a <- barplot(d,
legend.text=colnames(d),
col=color_palettes$response_palette[rownames(d)],
cex.names=namesCex,
cex.axis=axisCex,
cex.lab=labCex,
width=0.2,
xlim=c(0,1),
ylab="fraction of patients",
xaxt="n",
args.legend=list(x="topright",
bty="n",
fill=color_palettes$response_palette,
legend=c("CR/PR", "SD/PD"),
cex=legendCex))
text(x = a,
y = par("usr")[3] - 0.06,
labels = ifelse(colnames(d) == FALSE, "non-mutant", "mutant"),
srt = -45,
xpd = TRUE,
adj=0,
cex=namesCex,
xpd=TRUE)
mtext(nSamples,
side=3,
at=a,
line=0,
cex=0.7)
mtext("CCR w/o TP53, response",
side=3,
at=mean(a),
line=1,
font=titleF,
cex=titleCex)
par(mar=oldMar)
orgMar <- par()$mar
keggRes <- read.csv(file.path(tablesDir,
"SupplementaryTableS6.csv"),
stringsAsFactor=FALSE,
check.names=FALSE)
keggRes <- keggRes[, 1:6]
keggRes <- keggRes[!is.na(keggRes$S), ]
ind <- rev(1:nrow(keggRes))
keggRes <- keggRes[ind,]
up <- keggRes$Direction == "Up"
keggRes$"P (adj.)"[up] <- log10(keggRes$"P (adj.)"[up]) * -1
down <- keggRes$Direction == "Down"
keggRes$"P (adj.)"[down] <- log10(keggRes$"P (adj.)"[down])
resC1 <- read.csv(file.path(tablesDir,
"SupplementaryTableS3.csv"),
stringsAsFactor=FALSE,
check.names=FALSE)
keggRes$RepGenes <- NA_character_
for (i in 1:nrow(keggRes)) {
intGenes <- unlist(strsplit(keggRes$"Identified Genes"[i], ", "))
if (length(intGenes) > 7) {
intGenes <- resC1$Symbol[resC1$Symbol %in% intGenes][1:7]
intGenes <- c(intGenes, "...")
}
intGenes <- paste(intGenes, collapse=", ")
keggRes$RepGenes[i] <- intGenes
}
par(lheight=.7)
layout(matrix(c(1,0,0,0,1,0,0,0), 2, 4, byrow=TRUE))
proliferationPathways <- c("DNA replication",
"Cell cycle",
"Systemic lupus erythematosus",
"Alcoholism",
"Oocyte meiosis",
"Viral carcinogenesis",
"Progesterone-mediated oocyte maturation")
setCols <- ifelse(keggRes$Name %in% proliferationPathways,
"turquoise4",
ifelse(keggRes$Name %in% c("Fanconi anemia pathway", "Homologous recombination", "Mismatch repair", "Nucleotide excision repair", "Base excision repair", "p53 signaling pathway"),
"magenta4",
ifelse(keggRes$Name %in% "Cytokine-cytokine receptor interaction", "orange2", "darkgrey")))
a <- barplot(keggRes$"P (adj.)",
horiz=TRUE,
xlab="",
xlim=c(-1, 17),
col=setCols,
axes=FALSE)
mtext("KEGG pathway",
side=4,
at=max(a[,1]) + 1.5,
las=2,
line=3,
cex=titleCex)
mtext("Top-scoring genes",
side=4,
at=max(a[,1]) + 1.5,
las=2,
line=23,
cex=titleCex)
mtext(side=4,
at=a[,1],
text=paste0(keggRes$Name,
ifelse(keggRes$rep, "*", "")),
las=2,
line=1.5,
cex=(namesCex - 0.1),
col="black")
mtext(keggRes$RepGenes,
side=4,
at=a[,1],
line=22,
cex=(namesCex - 0.1),
las=2)
axis(side=1,
at=c(-2, seq(from=2, to=16, by=2)),
labels=c(2, seq(from=2, to=16, by=2)),
tick=TRUE,
xpd=TRUE,
cex=axisCex)
mtext(c("Up in responders", "Down"),
side=2,
line=1.5,
at=c(a[,1][10], a[,1][1]),
cex=axisCex)
segments(x0=-2.1,
x1=-2.1,
y0=a[,1][2] - 0.3,
y1=a[,1][length(a[,1])] + 0.3,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=-2.1,
x1=-2.1,
y0=a[,1][1] - 0.5,
y1=a[,1][1] + 0.5,
col="black",
xpd=TRUE,
lwd=1)
mtext(expression(paste("-log"["10"], " adj. p")),
side=1,
line=2.75,
at=7,
cex=axisCex)
legend("bottomright",
fill=c("turquoise4", "magenta4", "orange2"),
legend=c("Proliferation", "DNA damage response", expression(paste("TGF-", beta))),
bty="n",
cex=legendCex + 0.3)
segments(x0=0,
x1=0,
y0=-2,
y1=max(a[,1])+1,
col="black",
xpd=TRUE,
lwd=1)
par(mar=orgMar)
plotDots <- FALSE
orgMar <- par()$mar
tmpDat <- pData(cds2)
tmpDat$IFNG <- scale(m2["IFNG", ], center=TRUE, scale=TRUE)
tmpDat$IFNGR1 <- scale(m2["IFNGR1", ], center=TRUE, scale=TRUE)
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
for (sig in c("IFNG", "IFNGR1")) {
yLim1 <- -3
yLim2 <- ifelse(sig == "IFNG", 3, 5)
yLim <- c(yLim1, yLim2)
par(mar=c(5.5, 4.1, 2, 0.5))
pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P", sig, "score by binary Response:", pval))
yLab <- paste(sig, "expression")
a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
ylab=yLab,
#xlab="response",
col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
ylim=yLim,
whisklty = 1)
if (plotDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext(paste0(sig, ", response"),
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
x1=2,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3,
x1=4,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=1.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3.5,
x1=3.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=3.5,
y0=par("usr")[4] - yunit * 2,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
text(y=par("usr")[4] - yunit ,
x=2.5,
labels=pLevel(pval),
cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex))
}
## [1] "Wilcoxon test P IFNG score by binary Response: 5.2e-05"
## [1] "T test P IFNG score by binary Response: 9.1e-05"
## [1] "Wilcoxon test P IFNGR1 score by binary Response: 8.7e-05"
## [1] "T test P IFNGR1 score by binary Response: 0.00012"
par(mar=orgMar)
plotDots <- FALSE
binary <- FALSE
orgMar <- par()$mar
tmpDat <- pData(cds2)
tmpDat$TGFBR2 <- scale(m2["TGFBR2", ], center=TRUE, scale=TRUE)
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
sig <- "TGFBR2"
par(mar=c(5.5, 4.1, 2, 0.5))
pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))
## [1] "Wilcoxon test P TGFBR2 score by binary Response: 0.00029"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value,
2)
print(paste("T test P", sig, "score by binary Response:", pval))
## [1] "T test P TGFBR2 score by binary Response: 0.00019"
yLab <- paste(sig, "expression")
a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
ylab=yLab,
col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
cex.axis=axisCex,
cex.names=namesCex,
cex.lab=labCex,
ylim=c(-5, 5),
whisklty = 1)
if (plotDots) {
points(y=tmpDat[, sig],
x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
col=alpha("darkgrey", 0.7),
pch=16)
}
mtext(a$n,
side=3,
at=1:nlevels(tmpDat[, irf]),
line=0,
cex=0.75)
mtext(paste0(sig, ", response"),
side=3,
at=2.5,
line=1,
font=titleF,
cex=titleCex)
yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
x1=2,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3,
x1=4,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 4,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=1.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=3.5,
x1=3.5,
y0=par("usr")[4] - yunit * 4,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
segments(x0=1.5,
x1=3.5,
y0=par("usr")[4] - yunit * 2,
y1=par("usr")[4] - yunit * 2,
col="black",
xpd=TRUE,
lwd=1)
text(y=par("usr")[4] - yunit ,
x=2.5,
labels=pLevel(pval),
cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex))
par(mar=orgMar)
bio <- "TGFBR2"
tmpDat <- pData(cds2)
tmpDat$TGFBR2 <- m2["TGFBR2", ]
tmpDat <- tmpDat[!is.na(tmpDat[, bio]),]
tmpDat$group <- cut(tmpDat[, bio],
quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
labels=c("1st Quart",
"2nd Quart",
"3rd Quart",
"4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,
data=tmpDat)
plotSurvival2(survFit=fittedSurv,
survDiff=diffSurv,
diff.factor=as.factor(tmpDat$group),
main="TGFBR2, survival",
cols=rev(color_palettes[["irf_palette"]][1:4]),
ltypes=rep(1, nlevels(tmpDat$group)),
pval="none",
#coxphData=coxphDat,
plotMedian=TRUE,
plot.nrisk=FALSE,
xLab="OS (months)",
lwd=3,
cexMedian=axisCex - 0.1,
cexLegend=legendCex,
cex.lab=labCex)
print(coxph(Surv(os, censOS) ~
as.integer(tmpDat$group),
data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group),
## data = tmpDat)
##
## coef exp(coef) se(coef) z p
## as.integer(tmpDat$group) 0.13451 1.14398 0.05866 2.293 0.0218
##
## Likelihood ratio test=5.29 on 1 df, p=0.02147
## n= 348, number of events= 232
goi <- c("CD 8 T effector",
"gene19")
resp <- "binaryResponse"
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, resp]), ]
tmpDat[, resp] <- droplevels(tmpDat[, resp])
for (sig in goi) {
tmpDat[, sig] <- scale(tmpDat[, sig], center=TRUE, scale=TRUE)
}
tmpDat <- tmpDat[!is.na(tmpDat[, ml]) & tmpDat[, ml] > 0, ]
tmpDat[, ml] <- log2(tmpDat[, ml])
tmpDat <- tmpDat[! is.na(tmpDat$"Immune phenotype"),]
tmpDat$isExcluded <- as.factor(ifelse(tmpDat[, "Immune phenotype"] == "excluded",
"Yes", "No"))
varExp1 <- t(sapply(c(goi, ml), function(sig) {
fit <- glm(tmpDat[, resp] ~
tmpDat[, sig],
family="binomial")
fit$null.deviance - fit$deviance
}))
varExp1 <- setNames(varExp1[1,],
colnames(varExp1))
fit2.2 <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] + tmpDat[, "gene19"],
family="binomial")
varExp2.2 <- fit2.2$null.deviance - fit2.2$deviance
fit2.4 <- glm(tmpDat[, resp] ~
tmpDat[, ml] + tmpDat[, "CD 8 T effector"],
family="binomial")
varExp2.4 <- fit2.4$null.deviance - fit2.4$deviance
fit2.6 <- glm(tmpDat[, resp] ~
tmpDat[, ml] * tmpDat[, "gene19"],
family="binomial")
varExp2.6 <- fit2.6$null.deviance - fit2.6$deviance
varExp2 <- c(varExp2.2, varExp2.4, varExp2.6)
names(varExp2) <- c("Teff,TBRS", "Teff,TMB", "TBRS,TMB")
varExp2 <- sort(varExp2)
fit3.2 <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] + tmpDat[, ml] + tmpDat[, "gene19"] +
tmpDat[, ml]:tmpDat[, "gene19"] +
tmpDat[, "gene19"]*tmpDat$isExcluded,
family="binomial")
varExp3.2 <- fit3.2$null.deviance - fit3.2$deviance
names(varExp3.2) <- "Teff,TBRS,TMB"
allVars <- c(varExp1,
NA,
varExp2,
NA,
varExp3.2)
names(allVars) <- sub("CD 8 T effector", "Teff", names(allVars))
names(allVars) <- sub("gene19", "TBRS", names(allVars))
names(allVars) <- sub(ml, "TMB", names(allVars))
## repeat previous model to fit displayed data
fit_teff_0 <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"],
family="binomial")
fit_tgfb_0 <- glm(tmpDat[, resp] ~
tmpDat[, "gene19"],
family="binomial")
fit_teff_tgfb <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] +
tmpDat[, "gene19"],
family="binomial")
fit_tt_0 <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] +
tmpDat[, "gene19"],
family="binomial")
fit_ml_0 <- glm(tmpDat[, resp] ~
tmpDat[, ml],
family="binomial")
fit_tt_ml <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] +
tmpDat[, "gene19"] +
tmpDat[, ml],
family="binomial")
fit_tgfb_ml <- glm(tmpDat[, resp] ~
tmpDat[, "gene19"] +
tmpDat[, ml],
family="binomial")
fit_tgfb_ml_int <- glm(tmpDat[, resp] ~
tmpDat[, "gene19"] *
tmpDat[, ml],
family="binomial")
fit_teff_ml <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] +
tmpDat[, ml],
family="binomial")
fit_core_int <- glm(tmpDat[, resp] ~
tmpDat[, "CD 8 T effector"] +
tmpDat[, ml] +
tmpDat[, "gene19"] +
tmpDat[, "gene19"]:tmpDat[, ml] +
tmpDat[, "gene19"]*tmpDat$isExcluded,
family="binomial")
# likelihood ratio test Ps for display
tgfbP <- signif(anova(fit_tgfb_0, fit_teff_tgfb, test="Chisq")$"Pr(>Chi)"[2], 2)
teffP <- signif(anova(fit_teff_0, fit_teff_tgfb, test="Chisq")$"Pr(>Chi)"[2], 2)
ttP <- signif(anova(fit_tt_0, fit_tt_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
mlP <- signif(anova(fit_ml_0, fit_tt_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
teffmlP <- signif(anova(fit_ml_0, fit_teff_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
teffmlP2 <- signif(anova(fit_teff_0, fit_teff_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP <- signif(anova(fit_ml_0, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP2 <- signif(anova(fit_tgfb_0, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP3 <- signif(anova(fit_tgfb_ml, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP_int <- signif(anova(fit_tgfb_0, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP_int2 <- signif(anova(fit_ml_0, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)
tt_tgfb_mlP <- signif(anova(fit_tt_ml, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
core_interaction <- signif(anova(fit_tt_ml, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
core_tt <- signif(anova(fit_tt_0, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
core_tgfbml <- signif(anova(fit_tgfb_ml_int, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
core_teffml <- signif(anova(fit_teff_ml, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
pvals <- list(tgfbP=tgfbP,
teffP=teffP,
ttP=ttP,
mlP=mlP,
teffmlP=teffmlP,
teffmlP2=teffmlP2,
tgfbmlP=tgfbmlP,
tgfbmlP2=tgfbmlP2,
tgfbmlP3=tgfbmlP3,
tt_tgfb_mlP=tt_tgfb_mlP,
core_int=core_interaction,
tgfbmlP_int=tgfbmlP_int,
tgfbmlP_int2=tgfbmlP_int2,
core_tt=core_tt,
core_tgfbml=core_tgfbml,
core_teffml=core_teffml)
print(pvals)
## $tgfbP
## [1] 0.0032
##
## $teffP
## [1] 0.0026
##
## $ttP
## [1] 9e-07
##
## $mlP
## [1] 0.08
##
## $teffmlP
## [1] 0.2
##
## $teffmlP2
## [1] 4.9e-08
##
## $tgfbmlP
## [1] 0.23
##
## $tgfbmlP2
## [1] 6.6e-08
##
## $tgfbmlP3
## [1] 0.014
##
## $tt_tgfb_mlP
## [1] 0.057
##
## $core_int
## [1] 0.0052
##
## $tgfbmlP_int
## [1] 2.2e-08
##
## $tgfbmlP_int2
## [1] 0.024
##
## $core_tt
## [1] 1.9e-07
##
## $core_tgfbml
## [1] 0.016
##
## $core_teffml
## [1] 0.0028
pvals <- lapply(pvals, function(x) {
ifelse(x < 0.001, "***",
ifelse(x < 0.01, "**",
ifelse(x < 0.05, "*",
ifelse(x < 0.1, ".", "n.s."))))
})
oldMar <- par()$mar
par(mar=c(6, 3.9, 2, 0))
a <- barplot(allVars,
main="",
ylab="explained variance (%)",
col="grey",
width=0.06,
xlim=c(0,1),
ylim=c(0,58),
#args.legend=list(bg=alpha("white", 0.7),
# box.col=alpha("white", 0.7),
# cex=0.9),
xaxt="n",
cex.axis=0.9)
axis(1,
at=a[c(1:3, 5:7, 9)],
labels=FALSE)
text(x = a,
y=-3,
labels = names(allVars),
srt = -45,
xpd = TRUE,
adj=0,
cex=namesCex)
print("variance explained:")
## [1] "variance explained:"
print(allVars)
## Teff TBRS TMB Teff,TBRS Teff,TMB
## 4.125469 4.491943 32.246450 NA 13.169119 33.885674
## TBRS,TMB Teff,TBRS,TMB
## 39.742273 NA 50.042237
mtext(at = mean(a),
line=0.85,
side=3,
text="Variance explained, core pathways",
cex=titleCex,
font=titleF)
## improvement of Teff, TBRS over single models
text(x = a[5],
y=c(16.3, 14.8),
labels = unlist(pvals[c("teffP", "tgfbP")]),
#pos=3,
cex=0.9,
xpd=TRUE)
## improvement of Teff, TMB over single models
text(x = a[6],
y=c(37.3, 35.5),
labels = unlist(pvals[c("teffmlP2", "teffmlP")]),
#pos=3,
cex=0.9,
xpd=TRUE)
## improvement of TBRS, TMB over single models
text(x = a[7],
y=c(42.5, 41),
labels = unlist(pvals[c("tgfbmlP_int", "tgfbmlP_int2")]),
#pos=3,
cex=0.9,
xpd=TRUE)
## improvement of double on final core
segments(x0=a[7],
x1=a[9],
y0=51,
y1=51,
col="black",
xpd=TRUE,
lwd=2)
text(x = median(c(a[7], a[9])),
y=52.5,
labels = pvals[["core_tgfbml"]],
#pos=3,
cex=0.9,
xpd=TRUE)
segments(x0=a[6],
x1=a[9],
y0=55,
y1=55,
col="black",
xpd=TRUE,
lwd=2)
text(x = median(c(a[6], a[9])),
y=56.5,
labels = pvals[["core_teffml"]],
#pos=3,
cex=0.9,
xpd=TRUE)
segments(x0=a[5],
x1=a[9],
y0=59,
y1=59,
col="black",
xpd=TRUE,
lwd=2)
text(x = median(c(a[5], a[9])),
y=60.5,
labels = pvals[["core_tt"]],
#pos=3,
cex=0.9,
xpd=TRUE)
par(mar=oldMar)
[1] “2019-11-16 23:37:15” R version 3.5.3 (2019-03-11) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Mojave 10.14.4
Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages: [1] stats4 parallel grid stats graphics grDevices utils datasets [9] methods base
other attached packages: [1] corrplot_0.84 DT_0.9
[3] DESeq2_1.22.2 SummarizedExperiment_1.12.0
[5] DelayedArray_0.8.0 BiocParallel_1.16.6
[7] matrixStats_0.55.0 GenomicRanges_1.34.0
[9] GenomeInfoDb_1.18.2 IRanges_2.16.0
[11] S4Vectors_0.20.1 rmarkdown_1.16
[13] knitr_1.25 biomaRt_2.38.0
[15] readxl_1.3.1 edgeR_3.24.3
[17] limma_3.38.3 ROCR_1.0-7
[19] gplots_3.0.1.1 broom_0.5.2
[21] glmnet_2.0-18 foreach_1.4.7
[23] Matrix_1.2-17 MASS_7.3-51.4
[25] caret_6.0-84 lattice_0.20-38
[27] IMvigor210CoreBiologies_1.0.0 Biobase_2.42.0
[29] BiocGenerics_0.28.0 survminer_0.4.6
[31] survival_3.1-7 seriation_1.2-8
[33] viridis_0.5.1 viridisLite_0.3.0
[35] ggbeeswarm_0.6.0 ggpubr_0.2.3
[37] ggsci_2.9 GenomicDataCommons_1.6.0
[39] magrittr_1.5 PerformanceAnalytics_1.5.3
[41] xts_0.11-2 zoo_1.8-6
[43] ComplexHeatmap_1.20.0 forcats_0.4.0
[45] stringr_1.4.0 dplyr_0.8.3
[47] purrr_0.3.3 readr_1.3.1
[49] tidyr_1.0.0 tibble_2.1.3
[51] ggplot2_3.2.1 tidyverse_1.2.1
[53] workflowr_1.5.0
loaded via a namespace (and not attached): [1] utf8_1.1.4 tidyselect_0.2.5 htmlwidgets_1.5.1
[4] RSQLite_2.1.2 AnnotationDbi_1.44.0 TSP_1.1-7
[7] DESeq_1.34.1 munsell_0.5.0 codetools_0.2-16
[10] withr_2.1.2 colorspace_1.4-1 rstudioapi_0.10
[13] ggsignif_0.6.0 labeling_0.3 git2r_0.26.1
[16] GenomeInfoDbData_1.2.0 KMsurv_0.1-5 bit64_0.9-7
[19] rprojroot_1.3-2 vctrs_0.2.0 generics_0.0.2
[22] ipred_0.9-9 xfun_0.10 gclus_1.3.2
[25] R6_2.4.0 locfit_1.5-9.1 bitops_1.0-6
[28] assertthat_0.2.1 promises_1.1.0 scales_1.0.0
[31] nnet_7.3-12 beeswarm_0.2.3 gtable_0.3.0
[34] processx_3.4.1 timeDate_3043.102 rlang_0.4.1
[37] zeallot_0.1.0 genefilter_1.64.0 GlobalOptions_0.1.1
[40] splines_3.5.3 lazyeval_0.2.2 acepack_1.4.1
[43] ModelMetrics_1.2.2 checkmate_1.9.4 reshape2_1.4.3
[46] BiocManager_1.30.9 yaml_2.2.0 modelr_0.1.5
[49] crosstalk_1.0.0 backports_1.1.5 httpuv_1.5.2
[52] Hmisc_4.3-0 lava_1.6.6 tools_3.5.3
[55] ellipsis_0.3.0 RColorBrewer_1.1-2 Rcpp_1.0.3
[58] plyr_1.8.4 base64enc_0.1-3 progress_1.2.2
[61] zlibbioc_1.28.0 RCurl_1.95-4.12 ps_1.3.0
[64] prettyunits_1.0.2 rpart_4.1-15 GetoptLong_0.1.7
[67] cowplot_1.0.0 haven_2.2.0 cluster_2.1.0
[70] fs_1.3.1 data.table_1.12.6 circlize_0.4.8
[73] mime_0.7 hms_0.5.2 evaluate_0.14
[76] xtable_1.8-4 XML_3.98-1.20 gridExtra_2.3
[79] shape_1.4.4 compiler_3.5.3 KernSmooth_2.23-16
[82] crayon_1.3.4 htmltools_0.4.0 later_1.0.0
[85] Formula_1.2-3 geneplotter_1.60.0 lubridate_1.7.4
[88] DBI_1.0.0 rappdirs_0.3.1 cli_1.1.0
[91] quadprog_1.5-7 gdata_2.18.0 gower_0.2.1
[94] pkgconfig_2.0.3 km.ci_0.5-2 registry_0.5-1
[97] foreign_0.8-72 recipes_0.1.7 xml2_1.2.2
[100] annotate_1.60.1 vipor_0.4.5 XVector_0.22.0
[103] prodlim_2019.10.13 rvest_0.3.5 callr_3.3.2
[106] digest_0.6.22 cellranger_1.1.0 htmlTable_1.13.2
[109] survMisc_0.5.5 dendextend_1.12.0 curl_4.2
[112] shiny_1.4.0 gtools_3.8.1 rjson_0.2.20
[115] lifecycle_0.1.0 nlme_3.1-142 jsonlite_1.6
[118] fansi_0.4.0 pillar_1.4.2 fastmap_1.0.1
[121] httr_1.4.1 glue_1.3.1 iterators_1.0.12
[124] bit_1.1-14 class_7.3-15 stringi_1.4.3
[127] blob_1.2.0 latticeExtra_0.6-28 caTools_1.17.1.2
[130] memoise_1.1.0